In [1]:
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[1]:
In [2]:
import matplotlib.pyplot as plt
import numpy as np
import healpy as hp #1.40
import matplotlib.style
import matplotlib as mpl
mpl.style.use('classic')
from pylab import arange, show, cm
cmap = cm.bwr
cmap.set_under('w')
%matplotlib inline
params = {
'axes.labelsize': 24,
'font.size': 24,
'legend.fontsize': 20,
'xtick.labelsize': 22,
'ytick.labelsize': 22,
'lines.linewidth': 4,
'text.usetex': False,
# 'figure.autolayout': True,
'ytick.right': True,
'xtick.top': True,
'figure.figsize': [8, 6], # instead of 4.5, 4.5
'axes.linewidth': 2,
'xtick.major.size': 19,
'ytick.major.size': 19,
'xtick.minor.size': 10,
'ytick.minor.size': 10,
'xtick.major.width': 2,
'ytick.major.width': 2,
'xtick.minor.width': 2,
'ytick.minor.width': 2,
'xtick.major.pad': 10,
'ytick.major.pad': 10,
#'xtick.minor.pad': 14,
#'ytick.minor.pad': 14,
'xtick.direction': 'in',
'ytick.direction': 'in',
}
matplotlib.rcParams.update(params)
1. Input Data ¶
We will use Planck HFI 353 GHz map release 3: https://irsa.ipac.caltech.edu/data/Planck/release_3/all-sky-maps/
All sky maps are in HEALPix format, with Nside=2048 (for LFI 70GHz and HFI) with nested ordering.
COORDSYS= 'GALACTIC'
I, Q, U: The signal is given in units of Kcmb for 30-353 GHz.
II, QQ, UU, IQ, IU, QU: The signal in the unit of $K{cmb}^2$.
Beam: 4.818 arcmins
Unit Convention¶
From Kelvin to Jansky/sr
$ 1Jy = 10^{-26} W. m^{-2}. Hz ^{-1}$ ; $1W = kg. m^2.s^{−3} = \dfrac{1J}{s} = \dfrac{N.m}{s} = 1V.1A = ...$
sr = square radian
In the Rayleigh-Jeans limit (hν << kT), then
Where:
$𝑆_\nu$ is the flux density [$W m^{-2} Hz^{-1} sr^{-1}$],
c = 299792458 [m/s]
$\nu$ is the frequency [GHz],
${k_b}$ = 1.380649×10−23 [J⋅K−1] is Boltzmann's constant
$\theta_s$ is the size of the source (or the resolution in this case).
example:
beam_sigma = 50" # arcsec
beam_sigma = beam_sigma * $\pi$ / 60 /60 / 180 # radians
$\theta_s$ = beam_area = $2.\pi$.(beam_sigma)$^2$
As long as the source and wavelength are constant, brightness temperature and flux density can be interchanged.
ref0: https://arxiv.org/pdf/1111.1183.pdf
ref1: https://docs.astropy.org/en/stable/api/astropy.units.brightness_temperature.html
ref2: https://astronomy.stackexchange.com/questions/20704/how-are-people-converting-intensities-in-janskys-to-kelvin
In [3]:
kb = 1.380649*1e-23
c=299792458
nu=353*1e9
lamb = c**2 / nu**2 #m
K2Jy = 2*kb*nu**2/c**2 / 1e-26 # [J.K-1.m-2] = [W.s.K-1.m-2] = [W.Hz-1.m-2.K-1]
K2Jy
Out[3]:
In [4]:
planck_353 = hp.read_map("HFI_SkyMap_353_2048_R3.01_full.fits",field=(0,1,2),nest=False) #,3,4,5,6,7,8,9))
In [5]:
hp.nside2resol(2048, arcmin = True)
Out[5]:
In [6]:
I_stokes = planck_353[0]*K2Jy
Q_stokes = planck_353[1]*K2Jy
U_stokes = planck_353[2]*K2Jy
#hits = planck_353[3]
#II_cov = planck_353[4]*1e12 #uK^2
#IQ_cov = planck_353[5]*1e12
#IU_cov = planck_353[6]*1e12
#QQ_cov = planck_353[7]*1e12
#QU_cov = planck_353[8]*1e12
#UU_cov = planck_353[9]*1e12
In [7]:
plt.figure(figsize=(14,10)) ## maps in RING
hp.mollview(I_stokes,coord=["G"],unit='Jy',title="I Stokes",sub=[1,3,1],norm='hist',cmap=cmap)
hp.graticule()
hp.mollview(Q_stokes,coord=["G"],unit='Jy',title="Q Stokes",sub=[1,3,2],norm='hist',cmap=cmap)
hp.graticule()
hp.mollview(U_stokes,coord=["G"],unit='Jy',title="U Stokes",sub=[1,3,3],norm='hist',cmap=cmap)
hp.graticule()
2. Sky regions definition ¶
In [8]:
nside=2048
Ophiuchus_distance=140 #pc
l_ophi=354.0 # Galactic Longitude
b_ophi=17 # Galactic Latitude
delta_l_ophi,delta_b_ophi = 13.,13.
Lupus_distance=140. #pc
l_lupus= 340.0
b_lupus=12.7
delta_l_lupus,delta_b_lupus=12.,12.
Orion_distance = 450 #pc
l_orion = 212.0
b_orion = -16.0
delta_l_orion,delta_b_orion=16.,16.
l_taurus=172.5
b_taurus=-14.5
delta_l_taurus,delta_b_taurus=15.,15.
rot_ophi = [l_ophi,b_ophi]
rot_lupus = [l_lupus,b_lupus]
rot_orion = [l_orion,b_orion]
In [9]:
# uncomment these if you want to see locations
#hp.mollzoom(I_stokes,coord=["G"],title="I Stokes",rot=rot_ophi,norm='hist',cmap=cmap)
#hp.graticule()
#hp.mollzoom(I_stokes,coord=["G"],title="I Stokes",rot=rot_lupus,norm='hist',cmap=cmap)
#hp.graticule()
#hp.mollzoom(I_stokes,coord=["G"],title="I Stokes",rot=rot_orion,norm='hist',cmap=cmap)
#hp.graticule()
In [10]:
plt.figure(figsize=(14,10))
hp.gnomview(I_stokes,coord=["G"],title="I Stokes Ophiuchus",unit='Jy',rot=rot_ophi,sub=[1,3,1],xsize=4000,cmap=cmap,norm='hist')
hp.graticule()
hp.gnomview(I_stokes,coord=["G"],title="I Stokes Lupus",unit='Jy',rot=rot_lupus,sub=[1,3,2],xsize=4000,cmap=cmap,norm='hist')
hp.graticule()
hp.gnomview(I_stokes,coord=["G"],title="I Stokes Orion",unit='Jy',rot=rot_orion,sub=[1,3,3],xsize=2000,cmap=cmap,norm='hist')
hp.graticule()
3. Making masks ¶
In [11]:
def rectangle_mask_region(name,l,b,delta_l,delta_b,nside,write_mask=True,show_mask=True):
import healpy as hp
import numpy as np
# l : Galactic Longitude
# b : Galactic Latitude
theta = 90. - b # Spherical Latitude
phi = l # Spherical co-longitude
#calculate rectangle mask
theta_l,theta_u = theta - delta_b/2., theta + delta_b/2. # lower, upper
phi_l,phi_r = phi - delta_l/2., phi + delta_l/2. # left, right
v_ll = hp.ang2vec(np.radians(theta_l),np.radians(phi_l))
v_ul = hp.ang2vec(np.radians(theta_u),np.radians(phi_l))
v_ur = hp.ang2vec(np.radians(theta_u),np.radians(phi_r))
v_lr = hp.ang2vec(np.radians(theta_l),np.radians(phi_r))
vertices = np.vstack((v_ul, v_ll, v_lr, v_ur))
print('vertices',vertices)
indexes = hp.query_polygon(nside, vertices)
print('indexes',indexes)
t, p = np.degrees(hp.pix2ang(nside, indexes))
NPIX = hp.nside2npix(nside) ; print(NPIX)
m = np.arange(NPIX)*0. #m = np.zeros(hp.nside2npix(nside))
m[indexes]=1.
mask = m
if write_mask: hp.write_map('Mask'+name+str(nside)+'_.fits',mask,overwrite=True)
if show_mask:
hp.mollview(mask,coord=["G"], title="Mollview image RING")
hp.graticule()
return mask,indexes
In [12]:
mask_ophi,id_ophi=rectangle_mask_region('ophiuchus',l_ophi,b_ophi,delta_l_ophi,delta_b_ophi,nside,write_mask=True,show_mask=False)
mask_lupus,id_lupus=rectangle_mask_region('lupus',l_lupus,b_lupus,delta_l_lupus,delta_b_lupus,nside,write_mask=True,show_mask=False)
mask_orion,id_orion=rectangle_mask_region('orion',l_orion,b_orion,delta_l_orion,delta_b_orion,nside,write_mask=True,show_mask=False)
In [13]:
hp.gnomview(mask_lupus,coord=["G"],title="mask",rot=rot_lupus,sub=[1,3,3],xsize=600,norm='hist',cmap=cmap)
hp.graticule()
hp.mollview(mask_ophi+mask_lupus+mask_orion,coord=["G"],norm='hist',cmap=cmap)
hp.graticule()
The standard coordinates are the colatitude 𝜃, 0 at the North Pole, 𝜋/2 at the equator and 𝜋 at the South Pole and the longitude 𝜑 between 0 and 2𝜋 eastward, in a Mollview projection, 𝜑 = 0 is at the center and increases eastward toward the left of the map.
$\theta = \pi/2 - dec $
$\phi = ra $
$\theta = 90 - b$
$\phi = l$
4. Calculate selected regions ¶
In [14]:
I_ophiuchus = I_stokes*mask_ophi
Q_ophiuchus = Q_stokes*mask_ophi
U_ophiuchus = U_stokes*mask_ophi
I_orion = I_stokes*mask_orion
Q_orion = Q_stokes*mask_orion
U_orion = U_stokes*mask_orion
In [15]:
plt.figure(figsize=(14,10))
hp.gnomview(I_ophiuchus,title="I Ophiuchus",unit='Jy',rot=rot_ophi,sub=[1,3,1],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
hp.gnomview(Q_ophiuchus,title="Q Ophiuchus",unit='Jy',rot=rot_ophi,sub=[1,3,2],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
hp.gnomview(U_ophiuchus,title="U Ophiuchus",unit='Jy',rot=rot_ophi,sub=[1,3,3],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
In [16]:
#hp.mollzoom(I_ophiuchus,title="I Stokes",rot=rot_ophi,cmap=cmap,norm='hist')
#hp.graticule()
Polarization degrees:
Polarization angle is calculated:
In [17]:
p = np.sqrt(Q_ophiuchus**2 + U_ophiuchus**2) / I_ophiuchus
psi = 0.5*np.arctan(-U_ophiuchus/Q_ophiuchus) # confuse the sign: - or +
p_orion = np.sqrt(Q_orion**2 + U_orion**2) / I_orion
psi_orion = 0.5*np.arctan(-U_orion/Q_orion)
psi_orion
Out[17]:
In [18]:
# We can calculate just in pixels of Ophiuchus, it will avoid deviding by 0 warning!
#data = np.where(I_ophiuchus!=0)
#p = np.sqrt(Q_ophiuchus[id_ophi]**2 + U_ophiuchus[id_ophi]**2) / I_ophiuchus[id_ophi]
#psi = np.arctan(U_ophiuchus[id_ophi]/Q_ophiuchus[id_ophi])
In [19]:
plt.figure(figsize=(14,10))
hp.gnomview(p,coord=["G"],title="Polarization degrees",rot=rot_ophi,sub=[1,2,1],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
hp.gnomview(psi,coord=["G"],title="Polarization angle",rot=rot_ophi,sub=[1,2,2],xsize=480,norm='hist',cmap=cmap)
hp.graticule()
In [20]:
#m = np.arange(hp.nside2npix(nside))
#plt.plot(m, p)
#plt.xlabel('pixels')
#plt.ylabel('Polarization degrees')
fig, ax = plt.subplots(figsize=(10,7))
ax.hist(p,bins=40,histtype='step',linewidth=4.,color='b',range=[0,0.4],label='Ophiuchus')
ax.hist(p_orion,bins=40,histtype='step',linewidth=4.,color='g',range=[0,0.4],label='Orion')
ax.axvline(x=np.median(p[id_ophi]),linestyle='--',linewidth=4, color='b',label='Median:'+str(round(np.median(p[id_ophi]),3)))
ax.axvline(x=np.mean(p_orion[id_orion]),linestyle='-.',linewidth=4, color='g',label='Median:'+str(round(np.median(p_orion[id_orion]),3)))
plt.xlabel('Polarization degrees')
plt.ylabel('Occurance')
plt.legend()
plt.grid(linestyle='dotted')
fig, ax = plt.subplots(figsize=(10,7))
ax.hist(psi,bins=40,histtype='step',linewidth=4.,color='b',label='Ophiuchus')#,range=[0,0.6])
ax.hist(psi_orion,bins=40,histtype='step',linewidth=4.,color='g',label='Orion')#,range=[0,0.6])
ax.axvline(x=np.median(psi[id_ophi]),linestyle='--',linewidth=4, color='b',label='Median:'+str(round(np.median(psi[id_ophi]),3)))
ax.axvline(x=np.mean(psi_orion[id_orion]),linestyle='--',linewidth=4, color='g',label='Median:'+str(round(np.median(psi_orion[id_orion]),3)))
plt.xlabel('Polarization angle [rad]')
plt.ylabel('Occurance')
plt.legend(loc='lower center')
plt.grid(linestyle='dotted')
5. 2D real space data maps ¶
These are mollzoom plots to test the location we are going to take the data, they are just different in normalization color.
hp.mollzoom(I_stokes,coord=["G"],title="I Ophiuchus",unit='Jy',rot=rot_ophi,cmap=cmap,norm='hist')
hp.mollzoom(I_stokes,coord=["G"],title="I Ophiuchus",unit='Jy',rot=rot_ophi,cmap=cmap)
In [21]:
data = I_ophiuchus[id_ophi] # a linear dimension data
np.size(data)
np.size(id_ophi)
Out[21]:
In [22]:
histogram = plt.hist(data,bins=1000,range=[0,1e8])
In [23]:
NSIDE = 2048
hp.nside2resol(NSIDE, arcmin = True)
Out[23]:
In [24]:
N=np.size(data) #=197153 this is the number of pixels in a linear dimension
# since we are using lots of FFTs this should be a factor of 2^N
pix_size = hp.nside2resol(NSIDE, arcmin = True) #1.7177 # size of a pixel in arcminutes
planck_beam=4.8 #arc
## variables to set up the map plots; we have 197153 pixels, each pixel has 1.7177 arcminutes
X_width = N*pix_size/60. # horizontal map width in degrees , /60 = convert arcminutes to degrees
Y_width = N*pix_size/60. # vertical map width in degrees
#(we can devided by 2 for fast now)
print('X_width',X_width,'Y_width',Y_width)
In [25]:
import numpy as np
import matplotlib.pylab as plt
from pylab import arange, show, cm
cmap = cm.bwr
cmap.set_under('w')
from astropy.io import fits
hdulist = fits.open('HFI_SkyMap_353_2048_R3.01_full.fits')
hdulist.info()
#hdulist[1].header
In [26]:
## Test Interpolate values from healpy maps
#I_stokes = I_stokes*K2Jy # These maps in RING format, so if we use this we have to change hpx = HEALPix(nside=nside, order='RING', frame=Galactic())
#Q_stokes = I_stokes*K2Jy #
#U_stokes = I_stokes*K2Jy #
## Interpolate values from FITS file
I_stokes = hdulist[1].data['I_STOKES']*K2Jy # These data maps in NESTED format
Q_stokes = hdulist[1].data['Q_STOKES']*K2Jy #
U_stokes = hdulist[1].data['U_STOKES']*K2Jy #
In [27]:
# Set up the HEALPix projection of the data: the 353GHz map is in the Galactic coordinate
from astropy_healpix import HEALPix
from astropy.coordinates import Galactic
nside = hdulist[1].header['NSIDE'] # 2048
order = hdulist[1].header['ORDERING'] # NESTED
hpx = HEALPix(nside=nside, order=order, frame=Galactic()) # Make a map in Galactic frame ; NEST
5.1. In Galactic coordinate (longitude l,lattitude b) ¶
In [28]:
# Sample a X_width x Y_width grid which is coressponding to l,b, it depends on the data region.
from astropy import units as u
l=np.linspace( l_ophi - delta_l_ophi/2., l_ophi + delta_l_ophi/2., int(X_width)) * u.deg
b=np.linspace(b_ophi - delta_b_ophi/2., b_ophi + delta_l_ophi /2., int(Y_width)) * u.deg
l_grid,b_grid=np.meshgrid(l,b)
# Set up Astropy coordinate objects, the coordinate we want to see the objects
from astropy.coordinates import SkyCoord
coords = SkyCoord(l_grid.ravel(), b_grid.ravel(), frame='galactic')
In [29]:
Imap = hpx.interpolate_bilinear_skycoord(coords, I_stokes)
Qmap = hpx.interpolate_bilinear_skycoord(coords, Q_stokes)
Umap = hpx.interpolate_bilinear_skycoord(coords, U_stokes)
Imap_gal = np.reshape(Imap, ( int(X_width), int(Y_width) ) )[::-1,::-1] #swap x,y axis
Qmap_gal = np.reshape(Qmap, ( int(X_width), int(Y_width) ) )[::-1,::-1]
Umap_gal = np.reshape(Umap, ( int(X_width), int(Y_width) ) )[::-1,::-1]
In [30]:
#histogram = plt.hist(I_map,bins=1000,range=[0,1e8])
In [31]:
from matplotlib.colors import LogNorm
from mpl_toolkits.axes_grid1 import make_axes_locatable
ex=[l.value.min(), l.value.max(), b.value.min(), b.value.max() ]
plt.figure(figsize=(10,8) )
im = plt.imshow(Imap_gal, extent=ex , norm=LogNorm(), cmap=cmap) #,origin='left' )
plt.colorbar(im,pad=0.05)
plt.xlabel('Galactic longtitude $l$ [deg]')
plt.ylabel('Galactic latitude $b$ [deg]')
plt.show()
In [32]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
fig = plt.figure(figsize=(16, 12))
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=None)
ax1 = fig.add_subplot(121)
plt.xlabel(r"Galactic longtitude $l$ [deg]",fontsize="20")
plt.ylabel(r"Galactic latitude $b$ [deg]",fontsize="20")
im1 = ax1.imshow(Qmap_gal, extent=ex , norm=LogNorm(), cmap=cmap) #,origin='left' )
divider = make_axes_locatable(ax1)
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(im1, cax=cax, orientation='vertical')
ax2 = fig.add_subplot(122)
plt.xlabel(r"Galactic longtitude $l$ [deg]",fontsize="20")
plt.ylabel(r"Galactic latitude $b$ [deg]",fontsize="20")
im2 = ax2.imshow(Umap_gal,extent=ex, norm=LogNorm(), cmap=cmap)
divider = make_axes_locatable(ax2)
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(im2, cax=cax, orientation='vertical')
Out[32]:
In [33]:
outfile_I='I_stokes_Ophichus_Galactic_coord.fits' # HDU = Header Data Unit
outfile_Q='Q_stokes_Ophichus_Galactic_coord.fits' # HDU = Header Data Unit
outfile_U='U_stokes_Ophichus_Galactic_coord.fits' # HDU = Header Data Unit
hdu=fits.PrimaryHDU(Imap_gal)
hdu.writeto(outfile_I, clobber=True)
hdu=fits.PrimaryHDU(Qmap_gal)
hdu.writeto(outfile_Q, clobber=True)
hdu=fits.PrimaryHDU(Umap_gal)
hdu.writeto(outfile_U, clobber=True)
5.2. In equatorial coordinate (RA,DEC) ¶
In [34]:
from astropy import units as u
from astropy.coordinates import SkyCoord
c = SkyCoord(l_ophi, b_ophi, unit='deg', frame='galactic')
c.icrs
Out[34]:
In [35]:
ra_ophi = c.icrs.ra.hms
ra_ophi
Out[35]:
In [36]:
# Sample a X_width x Y_width grid which is coressponding to RA,DEC of region.
ra_ophi = c.icrs.ra.degree
dec_ophi = c.icrs.dec.degree
ra = np.linspace(ra_ophi - delta_l_ophi/2., ra_ophi + delta_l_ophi /2., int(X_width)) * u.deg
dec = np.linspace(dec_ophi - delta_b_ophi /2., dec_ophi + delta_l_ophi /2., int(Y_width)) * u.deg
ra_grid, dec_grid = np.meshgrid(ra, dec)
# Set up Astropy coordinate objects
from astropy.coordinates import SkyCoord
coords = SkyCoord(ra_grid.ravel(), dec_grid.ravel(), frame='icrs')
# Interpolate values
I_map = hpx.interpolate_bilinear_skycoord(coords, I_stokes)
Q_map = hpx.interpolate_bilinear_skycoord(coords, Q_stokes)
U_map = hpx.interpolate_bilinear_skycoord(coords, U_stokes)
Imap_equa = np.reshape(I_map, ( int(X_width), int(Y_width) ) )[::-1,::-1] #.swapaxes(-1,-1)#[...,::-1] #(-1,0)[...,::-2] #.swapaxes(-1,0)
Qmap_equa = np.reshape(Q_map, ( int(X_width), int(Y_width) ) )[::-1,::-1]
Umap_equa = np.reshape(U_map, ( int(X_width), int(Y_width) ) )[::-1,::-1]
In [37]:
from matplotlib.colors import LogNorm
ex=[ ra.value.max(),ra.value.min(),dec.value.min(), dec.value.max()]
plt.figure(figsize=(10,8))
im = plt.imshow(Imap_equa, extent=ex,norm=LogNorm(),cmap=cmap)#,origin='right')
plt.colorbar(im)
plt.xlabel('Right ascension [deg]')
plt.ylabel('Declination [deg]')
plt.show()
In [38]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
ex=[ ra.value.max(),ra.value.min(),dec.value.min(), dec.value.max()]
fig = plt.figure(figsize=(16, 12))
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=None)
ax1 = fig.add_subplot(121)
plt.xlabel('Right ascension [deg]')
plt.ylabel('Declination [deg]')
im1 = ax1.imshow(Qmap_equa, extent=ex , norm=LogNorm(), cmap=cmap) #,origin='left' )
divider = make_axes_locatable(ax1)
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(im1, cax=cax, orientation='vertical')
ax2 = fig.add_subplot(122)
plt.xlabel('Right ascension [deg]')
plt.ylabel('Declination [deg]')
im2 = ax2.imshow(Umap_equa,extent=ex, norm=LogNorm(), cmap=cmap)
divider = make_axes_locatable(ax2)
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(im2, cax=cax, orientation='vertical')
Out[38]:
In [39]:
outfile_I='I_stokes_Ophichus_RADEC_coord.fits' # HDU = Header Data Unit
outfile_Q='Q_stokes_Ophichus_RADEC_coord.fits' # HDU = Header Data Unit
outfile_U='U_stokes_Ophichus_RADEC_coord.fits' # HDU = Header Data Unit
hdu=fits.PrimaryHDU(Imap_equa)
hdu.writeto(outfile_I, clobber=True)
hdu=fits.PrimaryHDU(Qmap_equa)
hdu.writeto(outfile_Q, clobber=True)
hdu=fits.PrimaryHDU(Umap_equa)
hdu.writeto(outfile_U, clobber=True)
6. Test written fits data ¶
In [40]:
infile='I_stokes_Ophichus_Galactic_coord.fits'
In [41]:
hdu_list=fits.open(infile)
hdu_list.info()
In [42]:
I_ophi_data = hdu_list[0].data
print(type(I_ophi_data))
print(I_ophi_data.shape)
In [43]:
plt.imshow(I_ophi_data,norm=LogNorm(),cmap=cmap)
Out[43]:
In [ ]:
7. References ¶
[0] https://arxiv.org/pdf/1502.04123.pdf
[1] UNDERSTANDING POLARIZED DUST EMISSION FROM ρ OPHIUCHI A IN LIGHT OF GRAIN ALIGNMENT AND DISRUPTION BY RADIATIVE TORQUES https://arxiv.org/pdf/2007.10621.pdf
[2] Planck Early Results. XXV. Thermal dust in nearby molecular clouds https://arxiv.org/pdf/1101.2037.pdf
Appendix¶
Test coordinate transformation: https://en.wikipedia.org/wiki/Celestial_coordinate_system¶
In [44]:
# Set up Astropy coordinate objects
#from astropy.coordinates import SkyCoord
#coord = SkyCoord( frame='galactic', l="354.", b="17.",unit='deg')
In [45]:
from math import *
import numpy as np
def ga2equ(ga):
"""
Convert Galactic to Equatorial coordinates (J2000.0)
(use at own risk)
Input: [l,b] in decimal degrees
Returns: [ra,dec] in decimal degrees
Source:
- Book: "Practical astronomy with your calculator" (Peter Duffett-Smith)
- Wikipedia "Galactic coordinates"
Tests (examples given on the Wikipedia page):
>>> ga2equ([0.0, 0.0]).round(3)
array([ 266.405, -28.936])
>>> ga2equ([359.9443056, -0.0461944444]).round(3)
array([ 266.417, -29.008])
"""
l = radians(ga[0])
b = radians(ga[1])
# North galactic pole (J2000) -- according to Wikipedia
pole_ra = radians(192.859508)
pole_dec = radians(27.128336)
posangle = radians(122.932-90.0)
# North galactic pole (B1950)
#pole_ra = radians(192.25)
#pole_dec = radians(27.4)
#posangle = radians(123.0-90.0)
ra = atan2( (cos(b)*cos(l-posangle)), (sin(b)*cos(pole_dec) - cos(b)*sin(pole_dec)*sin(l-posangle)) ) + pole_ra
dec = asin( cos(b)*cos(pole_dec)*sin(l-posangle) + sin(b)*sin(pole_dec) )
return np.array([degrees(ra), degrees(dec)])
ra,dec=ga2equ([354.,17.])
ra,dec
Out[45]: